Sparse polynomial space approach to dissipative quantum systems: 
Application to the sub-ohmic spin-boson model 
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We propose a general numerical approach to open quantum systems with a coupling to bath 
degrees of freedom. The technique combines the methodology of polynomial expansions of spectral 
functions with the sparse grid concept from interpolation theory. Thereby we construct a Hilbert 
space of moderate dimension to represent the bath degrees of freedom, which allows us to perform 
highly accurate and efficient calculations of static, spectral and dynamic quantities using standard 
exact diagonalization algorithms. The strength of the approach is demonstrated for the phase 
transition, critical behaviour, and dissipative spin dynamics in the spin boson model. 
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Whenever a small quantum object, such as an atom, 
molecule or quantum dot, is not perfectly isolated it cou- 
ples to the degrees of freedom of its environment. In 
such an open quantum system the environment acts as a 
'bath' with which to exchange particles or energy with. 
A fermionic bath serves as a particle reservoir, while a 
bosonic bath accounts for dissipation Since the in- 
terest is only in the influence of the environment on the 
small quantum object, one may suspect that phenomeno- 
logical descriptions of open quantum systems, e.g. by 
Lindblad equations for dissipative baths, are sufficient. 
But in general correlations between the quantum system 
and the bath evolve, which can lead to strong renormal- 
ization as in the Kondo effect, or determine the time 
evolution of observables in unexpected ways. Simple phe- 
nomenological descriptions are obtained only within po- 
tentially unwarranted approximations such as weak cou- 
pling perturbation theory. To perform reliable compu- 
tations for open quantum systems including correlations 
with the environment is a challenging problem for theo- 
reticians. 

A generic and important example of an open quantum 
system is the spin-boson model [2|. Its Hamiltonian 



ibtb s 



(1) 



describes a spin- 1/2 (with Pauli matrices ai) coupled to 
a bosonic bath of oscillators, whose dynamics is given by 
Hb = T^i^ibt^i- The spin-boson coupling is specified 
by the spectral function 



J(u) 



\\ 2 l 8{uj-Lu l ) = -Lu 1 - s uj s <d{uj c -uj), (2) 



with a power-law dependence oc lu s up to a cutoff fre- 
quency uj c (we set ld c — 1 in the examples below). The 
spin-boson model shows rich physics beyond the dissi- 
pative spin dynamics at weak coupling. In the sub- 
ohmic (ohmic) regime s < 1 (s = 1) the model un- 
dergoes, for e = 0, a quantum phase transition (QPT) 



from a non-degenerate groundstate with zero magnetiza- 
tion m = (a z ) below a critical coupling a c = a c (A, s) 
to a two-fold degenerate groundstate with finite m ^ 
for a > a c . The existence of the QPT is a consequence 
of the coupling of the spin to bosons at low frequencies, 
which may entirely suppress the spin dynamics. In that 
respect the spin-boson model captures the renormaliza- 
tion aspect of Kondo physics. 

Only few methods are capable of accessing the QPT 
in the sub-ohmic spin-boson model. Among them we 
find powerful numerical techniques such as the Numer- 
ical Renormalization Group (NRG) 0, 0| or Quantum 
Monte Carlo (QMC) [|. Prominently missing in the 
above enumeration are techniques from the field of ex- 
act diagonalization (ED), which arc otherwise routinely 
used to yield highly accurate and unbiased results for 
strongly correlated systems Q. ED techniques require 
a finite-dimensional matrix representation of the model 
Hamiltonian. Once the matrix is given the Lanczos algo- 
rithm allows for the calculation of the groundstate and 
a few excited states, while Chebyshev expansion tech- 
niques such as the Kernel Polynomial Method (KPM) Q 
provide dynamic properties, e.g. spectral functions at 
zero or finite temperature, as well as the time-evolution 
of the wavefunction [8|. The main obstacle against this 
procedure for the spin-boson model and open quantum 
systems in general is that a finite Hamiltonian matrix 
involves discretization of the continuous spectral density 
J (to). Naive discretization, i.e. the approximate replace- 
ment of J(u>) by a sum of 5-peaks, requires cither a very 
large number of bosonic orbitals, which leads to matrices 
beyond any accessible size, or obtains results spoiled by 
discretization artefacts. 

The sparse polynomial space representation (SPSR) 
we propose in this Letter overcomes the ED restriction. 
It avoids the discretization of the bath spectral func- 
tion J(u>) and constructs a Hilbert space of moderate 
dimension to represent continuous bath degrees of free- 
dom with high resolution. In that way the SPSR extends 
the Chebyshev space method developed in Ref. , and it 
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FIG. 1: (Color online) Left panel: Two-dimensional sparse 
grid of level N g = 3 (circles) and N g = 4 (crosses). Right 
panel: Spectral function A(u>) = { j"; vac|5[w — H]\ f; vac) for 
A = 0, s = 0.5, a — 0.2 calculated using KPM. Keeping up 
to N b = 6 bosons, the SPSR to level N g = 10 contains 129284 
states. The comparable discrete grid (134596 states) contains 
only N p = 18 orbitals (here at equidistant energies Wj). 



becomes possible to perform efficient and accurate calcu- 
lations for open quantum systems using ED algorithms. 
As a non-trivial example we analyse the QPT and the 
dissipative spin dynamics in the spin-boson model. 

The Hilbert space of the Hamiltonian (Q} is the tensor 
product of the spin space C 2 with the bosonic Fock space 
B. To set up the SPSR for B we proceed in three steps: 
we (i) parametrize multiple bosonic excitations through 
symmetric wavefunctions as in first quantization, (ii) ex- 
pand these wavefunctions into orthogonal polynomials, 
(iii) select a sparse subspace of the polynomial space. 

For step (i) we fix an (unnormalized) density of states 
D{uS) = <5( w ~ w i) 011 [Oj w c ], which must be a smooth 
function for a continuous spectral function J(u>). In 
our numerics we use D(us) oc (1 — x 2 )" 1 / 2 with x = 
(2u>/uj c — 1) G (—1,1), which will lead to Chebyshev 
polynomials in step (ii). In first quantization any n- 
boson state \ip n ) is represented by a totally symmetric 
wavefunction ip n : [0,w c ]" — > C,Co >— > VVi(w). Here, the 
argument Co — {u>i, . . . , w n ) of the wavefunction gives the 
boson energies. We find that Hb multiplies the value 
ip n (Co) to argument Co by the total energy u>i. 

To express the Hamiltonian Eq. (JTJ) in our calcula- 
tions we further need the operators M + ) = Yli^i^i ■ 
These are bosonic operators up to normalization, since 
[b,b + ] — J2i^i — I dcoJ(uj). We choose the function 
\(lo) such that J{lo) = X(ui) 2 D(io), or X(u>i) = A, in 
comparison to Eq. ([1]). Then the single-boson state 
& + |vac) is represented by the wavefunction ip\{uj) = 
X{uS). Straightforward calculations show how to ob- 
tain the wavefunctions of any state b^ + '\ip n ). We note 
exemplarily, that for a single boson state |?/>i) with 
wavefunction the state b + \i/ji) has wavefunction 

-02(^1,^2) = (ipi(u;i)X(u>2) + X(u>i)ipi(u 2 ))/V2, while 
b\ipi) is the scalar J dLoD{u))\(<jj)ipi{to) . 

For step (ii), note that the scalar product of wavefunc- 



tions is given by 

(ipn,4>n) 



[0,« 



(3) 



Therefore we choose polynomials P m of degree m for m > 
subject to the orthonormality condition 



dujD{uj)Pi(uj)P m (uj) = 5i r< 



(4) 



For the above choice of D(u>), the P m are scaled and 
shifted Chebyshev polynomials. Any wavefunction ip n {Co) 
has an expansion 



i>n(&) = ^2 V'm Y[ P mi {^i) 



(5) 



in that complete polynomial function system. There- 
fore the multi-indices m enumerate the elements of an 
orthonormal basis of B. Instead with the wavefunction 
ip n (Co) we can calculate with the (totally symmetric) co- 
efficients V'm = /[ 0|U j„ dtiYliDfa^Pm^Ui)^^). 

Generally, ort hog onal polynomials P m obey a three- 
term recurrence [lfj] of the form P m +i = (a m co — b m )P m — 
c m P m -\. Owing to this recurrence the multiplication 
with J^i occurring for the operator Hb affects the co- 
efficients iprji only with index shifts by at most ±1. To 
obtain the operator M + ) we use the expansion X(u>) — 
Ylm^mPmi 1 ^) an d find, e.g., that 6 + |vac) has coeffi- 
cients ip m — X m . Similarly, for a single boson state 
\ip) with coefficients -0 m , the state b + \ip) has coefficients 

mi ,7712) = (ip mi ^m 2 + A mi -0m 2 )/v / 2, while b\ip) is the 
scalar J2 X m tp m . 

The bosonic Fock space and all relevant operators are 
now expressed by simple operations on a polynomial 
space. To prepare step (iii) notice that the selection of 
a finite dimensional subspace containing all polynomials 
up to degree N p is equivalent to naive discretization of 
J(lo), with Np + 1 energy levels Wi given as the zeroes of 
Pjv This discrete grid requires (n + N p )\/(n\N p \) 

coefficients to represent an n-boson state. To overcome 
the 'curse of dimension' expressed by the exponential 
growth of the binomial with n we resort to the concept 
of sparse grids ll| from interpolation theory. 



An n-dimensional sparse grid of level N g is a sub- 
set of the Cartesian grid with (2 Ng — 1)™ points (see 
Fig. [T]). With the sparse grid comes an interpolation 
formula that assigns a polynomial to given function val- 
ues at the sparse grid points. This interpolation has the 
property that functions of bounded variation are approxi- 
mated with high accuracy although the number of points 
is significantly smaller than in the Cartesian grid. For 
our purposes we do not access the points of the sparse 
grid directly. Instead we note that the sparse grid inter- 
polation formula is exact for a polynomial subspace of 
the full function space. Exactly this sparse polynomial 
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FIG. 2: (Color online) Left panels: Groundstate energy E as 
a function of oscillator shift 8; convergence of critical coupling 
a c with increasing Hilbert space size N g , number of bosons 
JVj,. Right panel: Phase diagram of the sub-ohmic spin-boson 
model for A = 0.1, i.e. a c as a function of s, in comparison 
to QMC/NRG data taken from Ref. d. 



space is selected for the SPSR. Assigning to a polynomial 
of degree m a logarithmic 'cost' co[m] — [\og 2 (m + 1)J 
(rounding down to an integer), we keep in step (iii) all 
polynomial basis states with multi-indices that satisfy 



co [mi] 

i=i 



< N n 



(6) 



For a single bosonic excitation (n — 1) the SPSR of level 
N g contains all polynomials with degree m < 2 Ng — 1. 
For n > 1, the SPSR contains only a small fraction of all 
polynomials, discarding those combinations where many 
polynomials have large degree. The motivation is that 
for multiple excitations the fine structure of the energy 
distributions among the various excitations becomes less 
important than for few excitations. Although the mo- 
tivation is related to Monte Carlo sampling of the state 
space, the SPSR is deterministic without statistical er- 
ror. Note further that, increasing N g , the SPSR is truly 
variational for the groundstate. 

In Fig. Q] the SPSR is compared to a discrete grid for 
the calculation of a spectral function. The discrete grid 
calculation is dominated by artefacts introduced by the 
inescapable restriction to a small number of orbitals. It 
is evident that the SPSR succeeds: Multiple bosonic ex- 
citations for continuous bath degrees of freedom are ac- 
curately represented with a moderate effort. Note that 
the SPRS resolves the jump discontinuity of A(u>) at the 
groundstate energy E s h = J dajJ(u>)/uj, and has uniform 
resolution over the full energy range. 

To put the SPSR to a severe test we calculate the 
phase transition in the sub-ohmic spin boson model. A 
NRG study of the QPT obtained for s < 1/2 critical 
behaviour incompatible with a mean-field transition ex- 
pected from the quantum-classical mapping to the Ising 
spin chain with long-range interactions [2j. Using QMC 
the authors recently corrected these findings Q , confirm- 
ing a mean-field transition. Apparently, the NRG calcu- 
lations of the critical behaviour suffered from a subtle 
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FIG. 3: (Color online) Upper left panel: Susceptibility x an d 
magnetization m as function of a, for s = 0.3 and A = 0.1. 
Upper right panel: Magnetization m as function of external 
field e, still for s = 0.3 and A = 0.1. The dashed straight 
lines indicate the slope of m for e — *• 0, which determines 
X- Lower panel: Critical behaviour of \ an d rn close to the 
phase transition, for A = 0.1 and s < 1/2. The curves for \ 
are multiplied with the indicated factors for better visibility. 
The straight lines indicate the critical behavior for a = (a — 
a c )/a c — > 0. The solid curves on the right show a fit to 
the ansatz m oc oP(l + O(q)), which results in the critical 
exponent j3 — 1/2 within numerical accuracy. 



error inherent to the renormalization scheme. In light of 
this controversy we use the SPSR to analyse the QPT 
independent of previous calculations. 

The QPT is best detected using the relation (bf+bi) = 
—2(Xi/uji)m between oscillator shift and magnetization 
in the groundstate. We therefore consider the Hamilto- 
nian 

H (5) = H + SJ2 *i (K + 2SE sh cr z + S 2 E sh , (7) 

i 

where the oscillator shift is introduced via the unitary 
transformation U(5) = exp[(5^ i (Ai/wi)(6^ — &<)]. In a 
certain sense U (S) prepares a classical mean-field state, 
while the quantum fluctuations are captured by the 
SPSR. Of course, the true groundstate energy E(6) of 
H(S) is independent of 5. But the SPSR becomes optimal 
if the oscillator shift, hence the average boson number, is 
small. Consequently, the numerical E(8) is minimal at fi- 
nite (zero) 5 if the true groundstate has finite (zero) mag- 
netization (see Fig. [2]). From E(6), calculated e.g. with 
the Lanczos algorithm, we obtain the critical coupling a c 
by simple bisection. Increasing the number of states in 
the SPSR the numerical values converge to the true a c 
(lower left panel) , which in turn yields the phase diagram 
(right panel) . In Fig. [3]we show the groundstate magneti- 
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FIG. 4: (Color online) Left panel: Magnetization m(t) = 
(i>(t)\a z \il>{t)) for A = 0.1, s = 1, a = 0.05 (with N g = 12, 
Nb = 6). For t < the system is prepared as a spin-up state 
with relaxed bosonic bath. The dashed curve shows the result 
from a discrete grid (N p — 25). Right panel: Decay of {cr x {t)) 
for A = 0, s = 1, a = 0.5. For t < the system is prepared 
as a spin singlet in the bosonic vacuum. We calculate a x (t) in 
the Heisenberg picture using a polynomial representation of 
operators. Already with N = 255 polynomials (corresponding 
to N g — 8) the numerical and analytical results match up to 
t = 1000. 

zation m and the susceptibility x — lim c _»o(<9m/de). The 
critical behaviour of the two quantities clearly confirms 
a mean-field transition for s < 0.5 with m ~ (a — aj 1 / 2 
and x ~ (o! c — a) -1 (Fig.|3J lower panel). Note that prob- 
ing for finite m or the divergence of x is an alternative 
to the above QPT criterion. The obtained values for a c 
agree with each other (cf. Fig. [2] and Fig. [3] for s = 0.3), 
but the above criterion is easier evaluated within the nu- 
merics, while e.g. x is obtained as a derivative. 

The analysis of the QPT demonstrates that the SPSR 
carries the unique virtues of ED techniques over to open 
quantum systems. Physical properties are found by the 
direct calculation of the corresponding observables. No 
scaling or extrapolation involving additional assumptions 
are required, no method specific quantities enter the dis- 
cussion. The computational effort is moderate, ranging 
from a few minutes to hours on standard PCs for the 
given results. Concerning their quality, our phase dia- 
gram is in perfect agreement with QMC and, taking the 
logarithmic NRG discretization into account, also with 
NRG. Our data for the critical behaviour confirm the 
QMC data, extrapolated to zero temperature. Here we 
can read off the critical behaviour directly from the nu- 
merical values. 

ED techniques have the overall advantage that, once 
the Hamiltonian matrix is given, almost any associated 
quantity can be obtained with high precision. Since the 
major interest is in the dynamics of open quantum sys- 
tem, we finally give a single example for the dissipative 
spin dynamics at weak coupling. Efficient time evolution 
with Chebyshev techniques [1, 0] gives the magnetization 
as a function of time (Fig.|U left panel). The curves are in 
perfec t ag reement with the results from time-dependent 
NRG [l2| and unitary perturbation theory [l]| . A spe- 
cial feature of our calculation is that it requires no ad- 



ditional damping, and no averaging over different bath 
discretizations. This results from the superior resolution 
provided by the SPSR even for multiple bosonic excita- 
tions. Although we have demonstrated that the SPSR is 
not restricted to weak coupling the time evolution close 
to the QPT deserves a careful examination that we post- 
pone to a future publication. To indicate the potential as 
the final example we show the decay of (a x (t)} for A = 0. 
For a finite number of polynomials the numerics exactly 
reproduces the analytical result, but only up to a finite 
time. With more polynomials that time can be easily 
made very large (Fig. [4] right panel). 

In conclusion, we introduced the SPSR as a novel ap- 
proach to static and dynamic properties of open quantum 
systems. The SPSR involves a highly accurate represen- 
tation of continuous bath degrees of freedom, which is 
based on the sparse grid concept applied to polynomial 
expansions of wavefunctions. It avoids the discretiza- 
tion artefacts that previously prevented the application 
of powerful ED techniques in presence of a bath. We 
demonstrated the strength of the SPSR for the QPT in 
the sub-ohmic spin-boson model, where we confirm the 
quantum-to-classical mapping for s < 1/2, and for the 
dissipative spin dynamics. Despite its current early state 
of development we believe to have presented the SPSR 
as a serious alternative to more established methods. 
An important issue for future work is the extension to 
fermionic baths, which is possible using antisymmetrized 
wavefunctions in step (i) of our construction. The effec- 
tiveness of SPRS in that case has yet to be assessed. 
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